Fonte Dati di Sopravvivenza

The Cancer Genome Atlas Program (TCGA) è un programma americano che ha caratterizzato 33 tipi di tumori e messo a disposizione della comunità scientifica i dati risultanti. L’obiettivo principale del progetto era creare un catalogo completo delle mutazioni genetiche e delle alterazioni molecolari responsabili dei vari tipi di cancro. Tra questi è incluso BRCA.

Linkedomics è un portale web che ha integrato i dati del TCGA e di un altro programma, e mette a disposizione dati già pre-processati.

BReast invasive CArcinoma

Requisiti per il tutorial

Per svolgere il seguente tutorial è necessario installare:

il software R disponibile per il proprio sistema operativo CRAN IT Mirror, e
l’IDE RStudio distribuito da posit.

Pacchetti R richiesti

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("TCGAbiolinks")
library("TCGAbiolinks")

require(reader)
require(dplyr)
require(janitor)
require(ggplot2)
require(survival)
require(survminer)
require(glmnet)

Download Dati Clinici

Il pacchetto R/Bioconductor TCGA Biolinks permette di scaricare direttamente in R i dataset presenti in Linkedomics.

# Query per i dati clinici
TCGA_BRCA_clinical = getLinkedOmicsData(
  project = "TCGA-BRCA",
  dataset = "Clinical"
)

dim(TCGA_BRCA_clinical)

## [1]   20 1098

Visualizzazione e Trasformazione

# Visualizziamo i nomi delle colonne (attributi)
TCGA_BRCA_clinical$attrib_name

##  [1] "years_to_birth"          "Tumor_purity"           
##  [3] "pathologic_stage"        "pathology_T_stage"      
##  [5] "pathology_N_stage"       "pathology_M_stage"      
##  [7] "histological_type"       "number_of_lymph_nodes"  
##  [9] "PAM50"                   "ER.Status"              
## [11] "PR.Status"               "HER2.Status"            
## [13] "gender"                  "radiation_therapy"      
## [15] "race"                    "ethnicity"              
## [17] "Median_overall_survival" "overall_survival"       
## [19] "status"                  "overallsurvival"

# Visualizziamo i dati di sopravvivenza (righe 18-20, per semplificare colonne 1-5)
TCGA_BRCA_clinical[18:20, 1:5]

## # A tibble: 3 × 5
##   attrib_name      TCGA.5L.AAT0 TCGA.5L.AAT1 TCGA.A1.A0SP TCGA.A2.A04V
##   <chr>            <chr>        <chr>        <chr>        <chr>       
## 1 overall_survival 1477         1471         584          1920        
## 2 status           0            0            0            1           
## 3 overallsurvival  1477,0       1471,0       584,0        1920,1

Nota. In caso di “overall_survival”, l’evento considerato nell’ analisi di sopravvivenza è la morte del paziente.

# Trasponiamo e puliamo i dati per avere un dataframe pazienti x dati clinici
TCGA_BRCA_clinical = t(TCGA_BRCA_clinical) %>%
  as.data.frame() %>%
  janitor::row_to_names(row_number = 1) %>%
  mutate(across(everything(), ~ type.convert(as.character(.), as.is = TRUE)))

dim(TCGA_BRCA_clinical)

## [1] 1097   20

Pulizia Dati

# Vediamo quanti NA abbiamo nei dati di sopravvvenza
print(paste("NA tempi sopravvivenza:", sum(is.na(TCGA_BRCA_clinical$overall_survival))))

## [1] "NA tempi sopravvivenza: 42"

print(paste("NA eventi:", sum(is.na(TCGA_BRCA_clinical$status))))

## [1] "NA eventi: 42"

# Rimuoviamo i pazienti con dati di sopravvivenza mancanti
TCGA_BRCA_clinical = TCGA_BRCA_clinical[
  !is.na(TCGA_BRCA_clinical$overall_survival) & 
    !is.na(TCGA_BRCA_clinical$status),]

dim(TCGA_BRCA_clinical)

## [1] 1055   20

Censoring plot

Il censoring plot ci mostra un summary dei dati di sopravvivenza, possiamo visualizzare i tempi (ascisse) e la presenza dell’evento o del censoring.

Usiamo solo i primi 30 pazienti per un plot più chiaro.

clin = TCGA_BRCA_clinical[1:30,]
clin = clin %>% arrange(overall_survival) # ordiniamo per tempi
ID = 1:nrow(clin)

clin %>%
  ggplot(
    aes(y = ID, x = overall_survival, shape = factor(status), color = factor(status))
  ) +
  geom_segment(aes(x = overall_survival, y = ID, xend = 0, yend = ID), linewidth = 2) +
  geom_point(size=3) +
  labs(x = "Days", y = "Patients") +
  scale_shape_discrete(name = "Status", labels = c("Censored","Dead")) +
  scale_color_manual(
    name = "Status",
    values = c("0" = "skyblue1", "1" = "red"),
    labels = c("0" = "Censored", "1" = "Dead")
  ) +
  scale_x_continuous(
    breaks = seq(0, max(clin$overall_survival, na.rm = TRUE), by = 180),    
    limits = c(0, max(clin$overall_survival, na.rm = TRUE))            
  ) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

Creazione Oggetto Surv

La funzione Surv() dal pacchetto survival combina i due vettori:

Time: Il tempo misurato durante lo studio clinico (di sopravvivenza o fino all’ultimo rilevamento).
Event: L’indicatore di stato (1 se l’evento è accaduto, 0 se il dato è censurato).

# Creare l'oggetto Surv
surv_object = Surv(
  time = TCGA_BRCA_clinical$overall_survival,
  event = TCGA_BRCA_clinical$status
)

# Visualizziamo i primi 20
head(surv_object, 20)

##  [1] 1477+ 1471+  584+ 1920  1099+ 2695+  595+ 3276+ 1165+  718+  954+  738+
## [13]  724+  661+  679+  343+  364+  336+ 1614+  883

L’output (es. 1477+) indica un paziente seguito per 1477 giorni e censurato (vivo all’ultimo follow-up). 883 indica un paziente deceduto (evento) al giorno 883.

Analisi Descrittiva

La curva di Kaplan-Meier è un grafico che mostra la probabilità di sopravvivenza stimata sull’asse Y in funzione del tempo sull’asse X.

Funzione a gradini: La curva è piatta tra un evento e l’altro e scende verticalmente solo nei momenti in cui si verificano uno o più eventi. L’altezza del “gradino” verso il basso dipende dal numero di eventi rispetto al numero di persone a rischio in quel momento.
Inizio al 100%: Al tempo 0, la probabilità di sopravvivenza è 1 (o 100%), poiché nessun evento si è ancora verificato.
Indicatori di censura: I tempi in cui i soggetti vengono censurati sono indicati sulla curva con piccoli segni per mostrare dove l’informazione è stata persa senza che si verificasse un evento.

Kaplan-Meier Globale

fit_km_global = survfit(
  surv_object ~ 1,
  data = TCGA_BRCA_clinical
)

ggsurvplot(
  fit_km_global,
  data = TCGA_BRCA_clinical,
  conf.int = TRUE,
  risk.table = TRUE,
  title = "Curva di Sopravvivenza Globale (TCGA-BRCA)",
  xlab = "Tempo (Giorni)"
)

Kaplan-Meier Stratificata

Creare curve separate per diversi gruppi ci permette di confrontare l’andamento della sopravvivenza rispetto una condizione (trattamento, stadiazione, sano vs malato).
Oltre ad una valutazione visiva, il log-rank test ci permette di stabilire se la differenza tra le curve è significativa.

La curva che rimane costantemente sopra l’ altra rappresenta il gruppo con la sopravvivenza migliore (cioè una minore probabilità di sperimentare l’evento in ogni dato momento).

Radioterapia

Stratifichiamo i pazienti rispetto il dato clinico radioterapia (yes/no).

#Rimuoviamo i pazienti con il dato clinico mancante
print(paste("NA radiation_therapy:", sum(is.na(TCGA_BRCA_clinical$radiation_therapy))))

## [1] "NA radiation_therapy: 79"

TCGA_BRCA_clinical_clean = TCGA_BRCA_clinical[!is.na(TCGA_BRCA_clinical$radiation_therapy), ]

fit_km_radiation_therapy = survfit(
  Surv(overall_survival, status) ~ radiation_therapy,
  data = TCGA_BRCA_clinical_clean
)

ggsurvplot(
  fit_km_radiation_therapy,
  conf.int = TRUE,
  data = TCGA_BRCA_clinical_clean,
  pval = TRUE,      #log-rank test
  risk.table = TRUE,   #tabella dei pazienti a rischio ad ogni intervallo di tempo
  legend.title = "Radioterapia"
)

Stadio patologico

print(paste("NA pathologic_stage:", sum(is.na(TCGA_BRCA_clinical$pathologic_stage))))

## [1] "NA pathologic_stage: 22"

TCGA_BRCA_clinical_clean2 = TCGA_BRCA_clinical[!is.na(TCGA_BRCA_clinical$pathologic_stage), ]

fit_km_path_stage = survfit(
  Surv(overall_survival, status) ~ pathologic_stage,
  data = TCGA_BRCA_clinical_clean2
)

ggsurvplot(
  fit_km_path_stage,
  conf.int = TRUE,
  data = TCGA_BRCA_clinical_clean2,
  risk.table = TRUE,
  legend.title = "Stadiazione"
)

Modelli di regressione: Cox Univariata

#Valutiamo l'impatto della variabile 'radiation_therapy'
fit_cox_rad = coxph(
  Surv(overall_survival, status) ~ radiation_therapy,
  data = TCGA_BRCA_clinical_clean
)

summary(fit_cox_rad)

## Call:
## coxph(formula = Surv(overall_survival, status) ~ radiation_therapy, 
##     data = TCGA_BRCA_clinical_clean)
## 
##   n= 976, number of events= 110 
## 
##                         coef exp(coef) se(coef)      z Pr(>|z|)   
## radiation_therapyyes -0.5519    0.5758   0.1914 -2.883  0.00394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                      exp(coef) exp(-coef) lower .95 upper .95
## radiation_therapyyes    0.5758      1.737    0.3957     0.838
## 
## Concordance= 0.57  (se = 0.028 )
## Likelihood ratio test= 8.32  on 1 df,   p=0.004
## Wald test            = 8.31  on 1 df,   p=0.004
## Score (logrank) test = 8.53  on 1 df,   p=0.004

Interpretazione: exp(coef) è l’Hazard Ratio (HR).
Un HR < 1 è protettivo verso l’evento. Un HR di 0.5758 significa che il gruppo “Radiation Therapy yes” ha un rischio di morte 0.5758 volte, in ogni istante, rispetto al gruppo “Radiation Therapy no”.

RNA sequencing

L’RNA-seq (RNA sequencing) è una tecnologia che permette di misurare l’espressione di tutti i geni rilevabili in una cellula, in un tessuto o in un organismo.

Come funziona (a grandi linee)?

Estrazione dell’RNA dalle cellule o dai tessuti (es. da campioni tumorali e/o sani).
Conversione in cDNA e sequenziamento (es. Illumina Hi-seq). Si ottengono milioni di “reads”.
Allineamento delle reads al genoma di riferimento.
1. Conteggio delle reads per gene -> si ottiene una matrice di espressione:
  
  Paziente 1 Paziente 2 Paziente 3 …
  
  gene 1 456 38 35
  
  gene 2 12 11 77
  
  …

	Paziente 1	Paziente 2	Paziente 3
gene 1	456	38	35
gene 2	12	11	77
…

Survival Analysis & RNA-seq

Dati RNA-seq: Misurano l’espressione di ~20.000 geni.
Obiettivo: Trovare geni la cui espressione (alta o bassa) è associata alla sopravvivenza.
Problema: Alta dimensionalità (p >> n).
Soluzione: Modelli di Cox penalizzati.

Download Dati RNA-seq

#Query per i dati di RNA-seq sequenziati con Illumina HiSeq
TCGA_BRCA_RNA = getLinkedOmicsData(
  project = "TCGA-BRCA",
  dataset = "RNAseq (HiSeq, Gene level)"
)

#Manipolazione del dataset per avere il dataframe campioni x geni
TCGA_BRCA_RNA = t(TCGA_BRCA_RNA) %>%
  as.data.frame() %>%
  janitor::row_to_names(row_number = 1) %>%
  mutate(across(everything(), ~ type.convert(as.character(.), as.is = TRUE)))

dim(TCGA_BRCA_RNA)

## [1]  1093 20155

# Selezioniamo solo i pazienti a cui sono associati dati di sopravvivenza
matching_ids = intersect(rownames(TCGA_BRCA_clinical), rownames(TCGA_BRCA_RNA))
BRCA_RNA = TCGA_BRCA_RNA[matching_ids,]
BRCA_clin = TCGA_BRCA_clinical[matching_ids,]

Modelli di regressione: Cox Multivariata

Gene signature

Genericamente prima di effettuare un’analisi multivariata si riducono i predittori con uno screening delle variabili in modo da avere \(d < p\) predittori informativi.

Selezioniamo alcuni geni che dalla letteratura sappiamo essere markers per il Breast Cancer: Breast Cancer gene 1 (BRCA1), Breast Cancer gene 2 (BRCA2), Estrogen Receptor Alpha (ESR1), Progesterone Receptor (PGR), Human Epidermal Growth Factor Receptor 2 (ERBB2)

gene_sign = c("BRCA1", "BRCA2", "ESR1", "PGR", "ERBB2")
BRCA_sign = BRCA_RNA[,gene_sign]

dim(BRCA_sign)

## [1] 1051    5

fit_cox_sign = coxph(
  Surv(BRCA_clin$overall_survival, BRCA_clin$status) ~ 
    BRCA1 + BRCA2 + ESR1 + PGR + ERBB2,
  data = BRCA_sign
)

summary(fit_cox_sign)

## Call:
## coxph(formula = Surv(BRCA_clin$overall_survival, BRCA_clin$status) ~ 
##     BRCA1 + BRCA2 + ESR1 + PGR + ERBB2, data = BRCA_sign)
## 
##   n= 1051, number of events= 151 
## 
##            coef exp(coef)  se(coef)      z Pr(>|z|)  
## BRCA1  0.120805  1.128404  0.086763  1.392   0.1638  
## BRCA2  0.019034  1.019216  0.077641  0.245   0.8063  
## ESR1   0.046387  1.047480  0.036998  1.254   0.2099  
## PGR   -0.069020  0.933308  0.033153 -2.082   0.0374 *
## ERBB2  0.006612  1.006634  0.054731  0.121   0.9038  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## BRCA1    1.1284     0.8862    0.9519     1.338
## BRCA2    1.0192     0.9811    0.8753     1.187
## ESR1     1.0475     0.9547    0.9742     1.126
## PGR      0.9333     1.0715    0.8746     0.996
## ERBB2    1.0066     0.9934    0.9042     1.121
## 
## Concordance= 0.592  (se = 0.028 )
## Likelihood ratio test= 8.27  on 5 df,   p=0.1
## Wald test            = 8.47  on 5 df,   p=0.1
## Score (logrank) test = 8.51  on 5 df,   p=0.1

Da questa analisi evinciamo che l’unico gene con p-value inferiore alla soglia convenzionale di \(\alpha=0.05\), è PGR. La sua espressione è statisticamente associata alla sopravvivenza globale.

In particolare, guardando a \(exp(coef)\), per ogni aumento unitario nell’espressione di PGR, il rischio di morte viene moltiplicato per 0.9333, mantenute costanti le altre variabili.

Indice Prognostico

L’indice prognostico è il predittore lineare PI = (coef_gene1 * expr_gene1) + (coef_gene2 * expr_gene2)…

merged_data = cbind(BRCA_sign, BRCA_clin)

prognostic_index = predict(
  fit_cox_sign,
  type = "lp", #linear prediction
  newdata = merged_data
)

# Stratifichiamo i pazienti in base alla mediana del PI
merged_data$risk_group = ifelse(
  prognostic_index > median(prognostic_index, na.rm=TRUE),
  "Alto Rischio",
  "Basso Rischio"
)

fit_km_risk = survfit(Surv(overall_survival, status) ~ risk_group, data = merged_data)

ggsurvplot(fit_km_risk,conf.int = TRUE, pval = TRUE, title = "Stratificazione per gruppi di rischio (Cox-Multivariata 5 Geni)")

Modelli di regressione: Cox Multivariata Penalizzata

glmnet

Perché glmnet?

Problema (p > n): Quando abbiamo più predittori (geni, \(p\)) che osservazioni (pazienti, \(n\)), la regressione di Cox standard fallisce.
Soluzione: Regressione Penalizzata. glmnet implementa vari tipi di penalizzazioni, tra cui:
- LASSO (alpha=1): Restringe i coefficienti e fa selezione delle variabili (molti coefficienti diventano zero).
- Ridge (alpha=0): Restringe i coefficienti ma mantiene tutte le variabili.
- Elastic Net (alpha tra 0 e 1): Un mix. Seleziona le variabili ed è più stabile di LASSO se i geni sono molto correlati.
\[ \hat{\beta} = \arg\max_{\beta} \left\{ \ell(\beta) - \lambda P(\beta) \right\} \]

\[ P(\beta) = \alpha \|\beta\|_1 + (1-\alpha)\|\beta\|_2^2 \]

Cox LASSO

Da uno screening delle variabili abbiamo identificato i geni più variabili (i.e. potenzialmente più informativi).

#Selezioniamo HV_genes (Highly Variable Genes)
HV_genes = read.csv("HV_genes.csv")$x
BRCA_HV = BRCA_RNA[, HV_genes]

# glmnet richiede una matrice per x e un oggetto Surv per y
x_matrix = as.matrix(BRCA_HV)
y_surv = Surv(BRCA_clin$overall_survival, BRCA_clin$status)

#glmnet applica LASSO impostando $alpha=1$
cvfit_lasso = cv.glmnet(
  x_matrix, y_surv,
  family = "cox", alpha = 1,
  nfolds = 5, # cv.glmnet trova il lambda ottimale con Cross-Validation (5-folds)
  type.measure = "C" # usiamo come metrica di valutazione il C-index
)

plot(cvfit_lasso)

#Variabili selezionate
coefficients_lasso = coef(cvfit_lasso, s="lambda.min")
selected_features = rownames(coefficients_lasso)[coefficients_lasso[,1] != 0]
length(selected_features)

## [1] 144

print(selected_features[1:20])

##  [1] "ACTL8"    "ADAM20"   "ADAMTS8"  "AKR1E2"   "ALS2CR12" "AMH"     
##  [7] "ASAH2"    "BAGE2"    "BCHE"     "BCL2L10"  "BCL8"     "C14orf50"
## [13] "C1orf168" "C1orf194" "C2orf39"  "C2orf82"  "C3orf45"  "C6orf142"
## [19] "C6orf81"  "C7orf54"

Cox Elastic Net

glmnet non ha un sistema interno per il tuning di \(alpha\) nel caso dell’ EN, quindi è possibile scrivere un ciclo per testare diversi alpha tra 0 e 1.

#alpha = seq(0.1, 1, length = 10)
#fitEN = list()
#for (i in 1:length(alpha)) {
#fitEN[[i]] = cv.glmnet(x_matrix, y_surv, family = "cox", alpha = alpha[i], nfolds = 5, type.measure = "C")
#}

# Troviamo il modello (alpha) con il C-index medio (cvm) migliore
#idx = which.max(sapply(fitEN, function(xx) {
#xx$cvm[xx$lambda == xx$lambda.min]
#}))
#cvfit_EN = fitEN[[idx]]
#best_alpha = alpha[idx]

Per semplificare impostiamo \(alpha=0.5\)

cvfit_EN = cv.glmnet(
  x_matrix, y_surv,
  family="cox", alpha=0.5,
  nfolds=5, type.measure="C"
)

plot(cvfit_EN)

coefficients_EN = coef(cvfit_EN, s="lambda.min")
selected_features_EN = rownames(coefficients_EN)[coefficients_EN[,1] != 0]

length(selected_features_EN)

## [1] 609

print(selected_features_EN[1:30])

##  [1] "ABCA4"    "ABCC6P1"  "ABCG4"    "ABP1"     "ACRV1"    "ACSBG2"  
##  [7] "ACTA1"    "ACTBL2"   "ACTC1"    "ACTL8"    "ACTN3"    "ADAM20"  
## [13] "ADAMTS19" "ADAMTS8"  "ADRA1A"   "AKNAD1"   "AKR1E2"   "ALDH3A1" 
## [19] "ALOXE3"   "ALS2CR11" "ALS2CR12" "AMAC1L2"  "AMH"      "AMN"     
## [25] "AMY1A"    "ANKRD2"   "ANO3"     "ANO5"     "ANXA13"   "APOF"

Confronto LASSO-Elastic Net

Calcoliamo i Prognostic Index (usando lambda.min, il lambda ottimale)

lp_lasso = predict(cvfit_lasso, newx=x_matrix, s="lambda.min", type="link")
lp_en = predict(cvfit_EN, newx=x_matrix, s="lambda.min", type="link")

# Stratifichiamo i gruppi di rischio
merged_data$risk_group_lasso = ifelse(lp_lasso > median(lp_lasso), "High", "Low")
merged_data$risk_group_enet  = ifelse(lp_en > median(lp_en), "High", "Low")

K-M (LASSO vs EN)

fit_km_lasso = survfit(Surv(overall_survival, status) ~ risk_group_lasso, data=merged_data)
ggsurvplot(fit_km_lasso, conf.int = TRUE, pval=TRUE, risk.table = TRUE, title="LASSO (alpha=1)", legend="none")

fit_km_enet = survfit(Surv(overall_survival, status) ~ risk_group_enet, data=merged_data)
ggsurvplot(fit_km_enet, conf.int = TRUE, pval=TRUE, risk.table = TRUE, title="Elastic Net (alpha=0.5)", legend="none")

Domande?

Grazie per l’attenzione!

Per approfondimenti e maggiori informazioni su metodi avanzati per l’analisi di sopravvivenza siete invitati alla Poster Session: Bridging Multi-omics Data ad Survival Analysis with Omics2Surv.

Esercizi

1 Stratifica i pazienti per classificazione istologica

Filtra i pazienti con dati mancanti per la variabile “PAM50”.
Fitta il modello.
Visualizza la kaplan meier stratificata.
Interpreta il log-rank test.

2 Costruisci un modello multivariato includendo oltre alla gene signature un dato clinico

Includi la variabile “radiation_therapy” alla cox multivariata.
Interpreta i risultati.
Stratifica per gruppi di rischio.

3 Effettua il tuning di alpha per Elastic Net

Prova a lanciare il codice per il tuning di alpha.

Ringraziamenti

This tutorial is part of the dissemination activities supported by the P2022BLN38 project Computational Approaches for the integration of multi-omics data – funded by European Union – Next Generation EU within the PRIN 2022 PNRR program (D.D. 1409 del 14-09-2022 Ministero dell’Università e della Ricerca) CUP B53D23027810001.

Analisi di Sopravvivenza in R

Un Caso Studio su Breast Cancer

Calanca Francesca, Raiconi Andrea, Angelini Claudia